#This project is statistical analysis for a 3X1 factorial experiment-1 and 3x2 factorial experiment-2. For experiment-1 contain four scenarios for each of three conditions.
# The conditions involve AI, Human and Human aided by AI. each of these are dummy coded # to run the regression model.
# TRUST, CREDEBILITY, AGREEABLENESS, AND FAIRNESS are the dependent variables for experiment-1. 
# TRUST, CREDEBILITY, QUALITY, AND BIAS are the dependent variables for experiment-2 as perceived by #participants. conditions is independent variable.

# Section-0: pre-process data to remove participants who does not qualify
# Section-1: regression analysis for experiment-1
# Section-2: MLM models experiment-1.
# Section-2: MLM models experiment-2.

########################################################################
#                                                                      #
#                     0 Data cleaning process                          #
#                                                                      #
########################################################################
#import relevant libraries
library(MASS)
library(dplyr)
library(scales)
library(performance)
library(reshape2)
library(nlme)
library(MuMIn)
library(car)
library(emmeans)
library(writexl)
library(haven)

library(writexl)
require(dplyr)
library(haven)

# loading the data from the first Poland study (the directory/data name goes below)
poland.data <- read_sav("xxxxxxxxxxxxxxxxxx") 

# loading the data from the first Poland study
poland.data <- read_sav("xxxxxxxxxxxxxxxxxx")

# filter participants who failed attention check and not active online
poland_final_data = poland.data %>% filter((filter1!=1 | filter2!=1 | filter3 !=1 | filter4!=1 ) & attn==2)
poland_final_data$survey_finish_time <- as.numeric(poland_final_data$survey_finish_time)
class(poland_final_data$survey_finish_time)

# Find 48% of the median time participants took to complete the experiment
median(as.numeric(poland.data$survey_finish_time))
quantile(poland.data$survey_finish_time)

#excluding-speeders i.e. participants who took less than 48% of the median time
#poland_final_data<-poland_final_data %>%filter(survey_finish_time > 337.92) #fist time data
poland_final_data<-poland_final_data %>%filter(survey_finish_time > 306)
nrow(poland_final_data)

#### excluding duplicates
temp <- poland_final_data %>% distinct(pid, .keep_all = T)
nrow(poland_final_data)


##################################################################################################
#                                                                                                #
#                   1 Regression analysis for experiment-1                                       #
#                                                                                                #
##################################################################################################
########################################################################
#                                                                      #
#                  1.1 Creating variables of interest                  #
#                                                                      #
########################################################################

#combining the three variables created for each question for three conditions for Scenario-1.

# Cond1FairS1, Cond2FairS1,Cond3FairS1 combined to make FairS1
poland_final_data$Cond1FairS1[is.na(poland_final_data$Cond1FairS1)]<-0 #replacing missing values with 0
poland_final_data$Cond2fairS1[is.na(poland_final_data$Cond2fairS1)]<-0
poland_final_data$cond3fairS1[is.na(poland_final_data$cond3fairS1)]<-0
poland_final_data$FAIR_S1 <- poland_final_data$Cond1FairS1 + poland_final_data$Cond2fairS1 + poland_final_data$cond3fairS1

# combining variables for BIAS
poland_final_data$Cond1BiasS1[is.na(poland_final_data$Cond1BiasS1)]<-0 #replacing missing values with 0
poland_final_data$Cond2biasS1[is.na(poland_final_data$Cond2biasS1)]<-0
poland_final_data$cond3biasS1[is.na(poland_final_data$cond3biasS1)]<-0
poland_final_data$BIAS_S1 <- poland_final_data$Cond1BiasS1 + poland_final_data$Cond2biasS1 + poland_final_data$cond3biasS1

# combining variables for TRUST
poland_final_data$Cond1trustS1[is.na(poland_final_data$Cond1trustS1)]<-0 #replacing missing values with 0
poland_final_data$cond2trustS1[is.na(poland_final_data$cond2trustS1)]<-0
poland_final_data$cond3trustS1[is.na(poland_final_data$cond3trustS1)]<-0
poland_final_data$TRUST_S1 <- poland_final_data$Cond1trustS1 + poland_final_data$cond2trustS1 + poland_final_data$cond3trustS1

# combining variables for LEGIT
poland_final_data$Cond1legitS1[is.na(poland_final_data$Cond1legitS1)]<-0 #replacing missing values with 0
poland_final_data$cond2legitS1[is.na(poland_final_data$cond2legitS1)]<-0
poland_final_data$cond3legitS1[is.na(poland_final_data$cond3legitS1)]<-0
poland_final_data$LEGIT_S1 <- poland_final_data$Cond1legitS1 + poland_final_data$cond2legitS1 + poland_final_data$cond3legitS1



# combining variabled for scenario-2
# Cond1FairS2, Cond2FairS2,Cond3FairS2 combined to make FairS2
poland_final_data$Cond1fairS2[is.na(poland_final_data$Cond1fairS2)]<-0 #replacing missing values with 0
poland_final_data$Cond2fairS2[is.na(poland_final_data$Cond2fairS2)]<-0
poland_final_data$Cond3fairS2[is.na(poland_final_data$Cond3fairS2)]<-0
poland_final_data$FAIR_S2 <- poland_final_data$Cond1fairS2 + poland_final_data$Cond2fairS2 + poland_final_data$Cond3fairS2

# combining variables for BIAS
poland_final_data$Cond1biasS2[is.na(poland_final_data$Cond1biasS2)]<-0 #replacing missing values with 0
poland_final_data$Cond2biasS2[is.na(poland_final_data$Cond2biasS2)]<-0
poland_final_data$Cond3BiasS2[is.na(poland_final_data$Cond3BiasS2)]<-0
poland_final_data$BIAS_S2 <- poland_final_data$Cond1biasS2 + poland_final_data$Cond2biasS2 + poland_final_data$Cond3BiasS2

# combining variables for TRUST
poland_final_data$Cond1TrustS2[is.na(poland_final_data$Cond1TrustS2)]<-0 #replacing missing values with 0
poland_final_data$Cond2trustS2[is.na(poland_final_data$Cond2trustS2)]<-0
poland_final_data$Cond3TrustS2[is.na(poland_final_data$Cond3TrustS2)]<-0
poland_final_data$TRUST_S2 <- poland_final_data$Cond1TrustS2 + poland_final_data$Cond2trustS2 + poland_final_data$Cond3TrustS2

# combining variables for LEGIT
poland_final_data$Cond1legitS2[is.na(poland_final_data$Cond1legitS2)]<-0 #replacing missing values with 0
poland_final_data$Cond2legitS2[is.na(poland_final_data$Cond2legitS2)]<-0
poland_final_data$Cond3legitS2[is.na(poland_final_data$Cond3legitS2)]<-0
poland_final_data$LEGIT_S2 <- poland_final_data$Cond1legitS2 + poland_final_data$Cond2legitS2 + poland_final_data$Cond3legitS2


# combining variabled for scenario-3

# Cond1FairS3, Cond2FairS3,Cond3FairS3 combined to make FairS3
poland_final_data$Cond1fairS3[is.na(poland_final_data$Cond1fairS3)]<-0 #replacing missing values with 0
poland_final_data$Cond2fairS3[is.na(poland_final_data$Cond2fairS3)]<-0
poland_final_data$Cond3fairS3[is.na(poland_final_data$Cond3fairS3)]<-0
poland_final_data$FAIR_S3 <- poland_final_data$Cond1fairS3 + poland_final_data$Cond2fairS3 + poland_final_data$Cond3fairS3

# combining variables for BIAS
poland_final_data$Cond1biasS3[is.na(poland_final_data$Cond1biasS3)]<-0 #replacing missing values with 0
poland_final_data$Cond2biasS3[is.na(poland_final_data$Cond2biasS3)]<-0
poland_final_data$Cond3biasS3[is.na(poland_final_data$Cond3biasS3)]<-0
poland_final_data$BIAS_S3 <- poland_final_data$Cond1biasS3 + poland_final_data$Cond2biasS3 + poland_final_data$Cond3biasS3

# combining variables for TRUST
poland_final_data$Cond1trustS3[is.na(poland_final_data$Cond1trustS3)]<-0 #replacing missing values with 0
poland_final_data$Cond2trustS3[is.na(poland_final_data$Cond2trustS3)]<-0
poland_final_data$Cond3trustS3[is.na(poland_final_data$Cond3trustS3)]<-0
poland_final_data$TRUST_S3 <- poland_final_data$Cond1trustS3 + poland_final_data$Cond2trustS3 + poland_final_data$Cond3trustS3

# combining variables for LEGIT
poland_final_data$Cond1legitS3[is.na(poland_final_data$Cond1legitS3)]<-0 #replacing missing values with 0
poland_final_data$Cond2legitS3[is.na(poland_final_data$Cond2legitS3)]<-0
poland_final_data$Cond3legitS3[is.na(poland_final_data$Cond3legitS3)]<-0
poland_final_data$LEGIT_S3 <- poland_final_data$Cond1legitS3 + poland_final_data$Cond2legitS3 + poland_final_data$Cond3legitS3

# combining variabled for scenario-4
# Cond1FairS4, Cond2FairS4,Cond3FairS4 combined to make FairS4
poland_final_data$Cond1fairS4[is.na(poland_final_data$Cond1fairS4)]<-0 #replacing missing values with 0
poland_final_data$Cond2fairS4[is.na(poland_final_data$Cond2fairS4)]<-0
poland_final_data$Cond3fairS4[is.na(poland_final_data$Cond3fairS4)]<-0
poland_final_data$FAIR_S4 <- poland_final_data$Cond1fairS4 + poland_final_data$Cond2fairS4 + poland_final_data$Cond3fairS4

# combining variables for BIAS
poland_final_data$Cond1biasS4[is.na(poland_final_data$Cond1biasS4)]<-0 #replacing missing values with 0
poland_final_data$Cond2biasS4[is.na(poland_final_data$Cond2biasS4)]<-0
poland_final_data$Cond3biasS4[is.na(poland_final_data$Cond3biasS4)]<-0
poland_final_data$BIAS_S4 <- poland_final_data$Cond1biasS4 + poland_final_data$Cond2biasS4 + poland_final_data$Cond3biasS4

# combining variables for TRUST
poland_final_data$Cond1trustS4[is.na(poland_final_data$Cond1trustS4)]<-0 #replacing missing values with 0
poland_final_data$Cond2trustS4[is.na(poland_final_data$Cond2trustS4)]<-0
poland_final_data$Cond3trustS4[is.na(poland_final_data$Cond3trustS4)]<-0
poland_final_data$TRUST_S4 <- poland_final_data$Cond1trustS4 + poland_final_data$Cond2trustS4 + poland_final_data$Cond3trustS4

# combining variables for LEGIT
poland_final_data$Cond1legitS4[is.na(poland_final_data$Cond1legitS4)]<-0 #replacing missing values with 0
poland_final_data$Cond2legitS4[is.na(poland_final_data$Cond2legitS4)]<-0
poland_final_data$Cond3legitS4[is.na(poland_final_data$Cond3legitS4)]<-0
poland_final_data$LEGIT_S4 <- poland_final_data$Cond1legitS4 + poland_final_data$Cond2legitS4 + poland_final_data$Cond3legitS4





#### creating dummy and mean variables for study1 data
# poland_final_data$AI.dummy <- ifelse(poland_final_data$CONDITION %in% c(2,3),0,1)
# poland_final_data$Human.dummy <- ifelse(poland_final_data$CONDITION %in% c(1,3),0,1)
# poland_final_data$AI_Human.dummy <- ifelse(poland_final_data$CONDITION %in% c(1,2),0,1)

#### creating dummy and mean variables for study2 data
poland_final_data$AI.dummy <- ifelse(poland_final_data$conditionexp1 %in% c(2,3),0,1)
poland_final_data$Human.dummy <- ifelse(poland_final_data$conditionexp1 %in% c(1,3),0,1)
poland_final_data$AI_Human.dummy <- ifelse(poland_final_data$conditionexp1 %in% c(1,2),0,1)

poland_final_data$pl1mean.trust <- (poland_final_data$TRUST_S1+poland_final_data$TRUST_S2+ poland_final_data$TRUST_S3+ poland_final_data$TRUST_S4)/4
poland_final_data$pl1mean.fair <- (poland_final_data$FAIR_S1+poland_final_data$FAIR_S2+ poland_final_data$FAIR_S3+ poland_final_data$FAIR_S4)/4
poland_final_data$pl1mean.biased <- (poland_final_data$BIAS_S1+poland_final_data$BIAS_S2+ poland_final_data$BIAS_S3+ poland_final_data$BIAS_S4)/4
poland_final_data$pl1mean.legit <- (poland_final_data$LEGIT_S1+poland_final_data$LEGIT_S2+ poland_final_data$LEGIT_S3+ poland_final_data$LEGIT_S4)/4


write_xlsx(poland_final_data, "poland_final_data.xlsx")


########################################################################
#                                                                      #
# 1.2  Randomization check across different control variables          #
#                                                                      #
########################################################################

### randomization check on Income, Gender, Education, Age and Ideology

# for study1 data
# summary(aov(scale(CONDITION)~scale(sex) , data = subset(poland_final_data,!is.na(sex))), na.rm=T)
# summary(aov(scale(CONDITION)~scale(education), data = subset(poland_final_data, !is.na(education))), na.rm=T)
# summary(aov(scale(CONDITION)~scale(age), data = subset(poland_final_data, !is.na(age))), na.rm=T)
# summary(aov(scale(CONDITION)~scale(ideology), data = subset(poland_final_data, !is.na(ideology))), na.rm=T)

# for study2 data
summary(aov(scale(conditionexp1)~scale(sex) , data = subset(poland_final_data,!is.na(sex))), na.rm=T)
summary(aov(scale(conditionexp1)~scale(edu), data = subset(poland_final_data, !is.na(edu))), na.rm=T)
summary(aov(scale(conditionexp1)~scale(age), data = subset(poland_final_data, !is.na(age))), na.rm=T)
summary(aov(scale(conditionexp1)~scale(ideology), data = subset(poland_final_data, !is.na(ideology))), na.rm=T)


########################################################################
#                                                                      #
# 1.3 Regression analysis for study1 (aggregated dependent variables)  #
#                                                                      #
########################################################################
summary(aov(scale(pl1mean.trust) ~ scale(conditionexp1), data = poland_final_data), na.rm=T)
summary(aov(scale(pl1mean.fair) ~ scale(conditionexp1) , data = poland_final_data), na.rm=T)
summary(aov(scale(pl1mean.biased) ~ scale(conditionexp1) , data = poland_final_data), na.rm=T)
summary(aov(scale(pl1mean.legit) ~ scale(conditionexp1), data = poland_final_data), na.rm=T)

summary(lm(scale(pl1mean.trust) ~ scale(AI.dummy) + scale(AI_Human.dummy) , data = poland_final_data), na.rm=T)
summary(lm(scale(pl1mean.fair) ~ scale(AI.dummy) + scale(AI_Human.dummy) , data = poland_final_data), na.rm=T)
summary(lm(scale(pl1mean.biased) ~ scale(AI.dummy) + scale(AI_Human.dummy) , data = poland_final_data), na.rm=T)
summary(lm(scale(pl1mean.legit) ~ scale(AI.dummy) + scale(AI_Human.dummy), data = poland_final_data), na.rm=T)

summary(lm(scale(pl1mean.trust) ~ scale(Human.dummy) + scale(AI_Human.dummy) , data = poland_final_data), na.rm=T)
summary(lm(scale(pl1mean.fair) ~ scale(Human.dummy) + scale(AI_Human.dummy) , data = poland_final_data), na.rm=T)
summary(lm(scale(pl1mean.biased) ~ scale(Human.dummy) + scale(AI_Human.dummy) , data = poland_final_data), na.rm=T)
summary(lm(scale(pl1mean.legit) ~ scale(Human.dummy) + scale(AI_Human.dummy), data = poland_final_data), na.rm=T)


########################################################################
#                                                                      #
#  1.4  Regression analysis for study1 (Scenariowise for each dv)      #
#                                                                      #
########################################################################
################################################### Scenario-1 ######################################################


summary(lm((TRUST_S1) ~ (AI.dummy) + (AI_Human.dummy) , data = poland_final_data), na.rm=T)
summary(lm((FAIR_S1) ~ (AI.dummy) + (AI_Human.dummy) , data = poland_final_data), na.rm=T)
summary(lm((BIAS_S1) ~ (AI.dummy) + (AI_Human.dummy) , data = poland_final_data), na.rm=T)
summary(lm((LEGIT_S1) ~ (AI.dummy) + (AI_Human.dummy), data = poland_final_data), na.rm=T)

summary(lm((TRUST_S1) ~ (Human.dummy) + (AI_Human.dummy) , data = poland_final_data), na.rm=T)
summary(lm((FAIR_S1) ~ (Human.dummy) + (AI_Human.dummy) , data = poland_final_data), na.rm=T)
summary(lm((BIAS_S1) ~ (Human.dummy) + (AI_Human.dummy) , data = poland_final_data), na.rm=T)
summary(lm((LEGIT_S1) ~ (Human.dummy) + (AI_Human.dummy), data = poland_final_data), na.rm=T)

################################################### Scenario-2 ######################################################

summary(lm((TRUST_S2) ~ (AI.dummy) + (AI_Human.dummy), data = poland_final_data), na.rm=T)
summary(lm((FAIR_S2) ~ (AI.dummy) + (AI_Human.dummy), data = poland_final_data), na.rm=T)
summary(lm((BIAS_S2) ~ (AI.dummy) + (AI_Human.dummy) , data = poland_final_data), na.rm=T)
summary(lm((LEGIT_S2) ~ (AI.dummy) + (AI_Human.dummy), data = poland_final_data), na.rm=T)

summary(lm((TRUST_S2) ~ (Human.dummy) + (AI_Human.dummy), data = poland_final_data), na.rm=T)
summary(lm((FAIR_S2) ~ (Human.dummy) + (AI_Human.dummy) , data =poland_final_data), na.rm=T)
summary(lm((BIAS_S2) ~ (Human.dummy) + (AI_Human.dummy), data = poland_final_data), na.rm=T)
summary(lm((LEGIT_S2) ~ (Human.dummy) + (AI_Human.dummy), data = poland_final_data), na.rm=T)

################################################### Scenario-3 ######################################################

summary(lm((TRUST_S3) ~ (AI.dummy) + (AI_Human.dummy) , data = poland_final_data), na.rm=T)
summary(lm((FAIR_S3) ~ (AI.dummy) + (AI_Human.dummy) , data = poland_final_data), na.rm=T)
summary(lm((BIAS_S3) ~ (AI.dummy) + (AI_Human.dummy) , data = poland_final_data), na.rm=T)
summary(lm((LEGIT_S3) ~ (AI.dummy) + (AI_Human.dummy), data = poland_final_data), na.rm=T)

summary(lm((TRUST_S3) ~ (Human.dummy) + (AI_Human.dummy) , data = poland_final_data), na.rm=T)
summary(lm((FAIR_S3) ~ (Human.dummy) + (AI_Human.dummy) , data = poland_final_data), na.rm=T)
summary(lm((BIAS_S3) ~ (Human.dummy) + (AI_Human.dummy) , data = poland_final_data), na.rm=T)
summary(lm((LEGIT_S3) ~ (Human.dummy) + (AI_Human.dummy), data = poland_final_data), na.rm=T)

################################################### Scenario-4 ######################################################

summary(lm((TRUST_S4) ~ (AI.dummy) + (AI_Human.dummy) , data = poland_final_data), na.rm=T)
summary(lm((FAIR_S4) ~ (AI.dummy) + (AI_Human.dummy) , data = poland_final_data), na.rm=T)
summary(lm((BIAS_S4) ~ (AI.dummy) + (AI_Human.dummy) , data = poland_final_data), na.rm=T)
summary(lm((LEGIT_S4) ~ (AI.dummy) + (AI_Human.dummy), data = poland_final_data), na.rm=T)

summary(lm((TRUST_S4) ~ (Human.dummy) + (AI_Human.dummy) , data = poland_final_data), na.rm=T)
summary(lm((FAIR_S4) ~ (Human.dummy) + (AI_Human.dummy) , data = poland_final_data), na.rm=T)
summary(lm((BIAS_S4) ~ (Human.dummy) + (AI_Human.dummy) , data = poland_final_data), na.rm=T)
summary(lm((LEGIT_S4) ~ (Human.dummy) + (AI_Human.dummy), data = poland_final_data), na.rm=T)

########################################################################
#                                                                      #
#         1.6 Internal consistency check and factor analysis            #
#                             Appendix-5 and 7                         #
########################################################################
# calculating Cronbach's alpha for internal consistency for PL sample
# use the subset of the questions whose internal consistency is to be calculated
library(psych)

# for study 2
#alpha_s1s1 <- poland_final_data[, 127:130]
# alpha(alpha_s1s1)
# 
# alpha_s1s2 <- poland_final_data[, 131:134]
# alpha(alpha_s1s2)
# 
# alpha_s1s3 <- poland_final_data[, 135:138]
# alpha(alpha_s1s3)
# 
# alpha_s1s4 <- poland_final_data[, 139:142]
# alpha(alpha_s1s4)

#for study 1
alpha_s1s1 <- poland_final_data[, 84:87]
alpha(alpha_s1s1)

alpha_s1s2 <- poland_final_data[, 88:91]
alpha(alpha_s1s2)

alpha_s1s3 <- poland_final_data[, 92:95]
alpha(alpha_s1s3)

alpha_s1s4 <- poland_final_data[, 96:99]
alpha(alpha_s1s4)


# 2.Choose the "right" number of factors to retain
##parallel analysis
library(nFactors)
library(psych)
library(GPArotation)

correl = cor(alpha_s1s1, use = "pairwise.complete.obs")
correl
symnum(correl)

ev <- eigen(cor(alpha_s1s1)) # get eigenvalues 
ap <- parallel(subject=nrow(x = alpha_s1s1),var=ncol(x = alpha_s1s1), rep=100, cent=.05) 
nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea) 
plotnScree(nS)
ap$eigen
nofactors1 <- fa.parallel(alpha_s1s1, fm='ml', fa= 'fa')
sum(nofactors1$fa.values >.7)


s1.fa1 <- fa(r = cor(alpha_s1s1), nfactors = 2, rotate = "oblimin", fm = "ml")
plot(s1.fa1,labels=names(x = alpha_s1s1),cex=.7, ylim=c(-.1,1)) 
s1.fa1

############### scenario2############
correl = cor(alpha_s1s2, use = "pairwise.complete.obs")
symnum(correl)

# 2.Choose the "right" number of factors to retain
##parallel analysis

ev <- eigen(cor(alpha_s1s2)) # get eigenvalues 
ap <- parallel(subject=nrow(x = alpha_s1s2),var=ncol(x = alpha_s1s2), rep=100, cent=.05) 
nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea) 
plotnScree(nS)
ap$eigen

nofactors2 <- fa.parallel(alpha_s1s2, fm='ml', fa= 'fa')
sum(nofactors2$fa.values >.7)


s1.fa2 <- fa(r = cor(alpha_s1s2), nfactors = 2, rotate = "oblimin", fm = "ml")
plot(s1.fa2,labels=names(x = alpha_s1s2),cex=.7, ylim=c(-.1,1)) 
s1.fa2

############### scenario3 ############
correl = cor(alpha_s1s3, use = "pairwise.complete.obs")
symnum(correl)

# 2.Choose the "right" number of factors to retain
##parallel analysis
library(nFactors) 
ev <- eigen(cor(alpha_s1s3)) # get eigenvalues 
ap <- parallel(subject=nrow(x = alpha_s1s3),var=ncol(x = alpha_s1s3), rep=100, cent=.05) 
nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea) 
plotnScree(nS)
ap$eigen

nofactors3 <- fa.parallel(alpha_s1s3, fm='ml', fa= 'fa')
sum(nofactors3$fa.values >.7)


s1.fa3 <- fa(r = cor(alpha_s1s3), nfactors = 2, rotate = "oblimin", fm = "ml")
plot(s1.fa3,labels=names(x = alpha_s1s3),cex=.7, ylim=c(-.1,1)) 
s1.fa3


########################## Scenario4 #######################

correl = cor(alpha_s1s4, use = "pairwise.complete.obs")
symnum(correl)

# 2.Choose the "right" number of factors to retain
##parallel analysis

ev <- eigen(cor(alpha_s1s4))# get eigenvalues 
ev$values
ap <- parallel(subject=nrow(x = alpha_s1s4),var=ncol(x = alpha_s1s4), rep=100, cent=.05) 
nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea) 
plotnScree(nS)
ap$eigen

nofactors4 <- fa.parallel(alpha_s1s4, fm='ml', fa= 'fa')
sum(nofactors4$fa.values >.7)


s1.fa4 <- fa(r = cor(alpha_s1s4), nfactors = 2, rotate = "oblimin", fm = "ml")
plot(s1.fa4,labels=names(x = alpha_s1s4),cex=.7, ylim=c(-.1,1)) 
s1.fa4



##################################################################################################
#                                                                                                #
#                   2  Multilevel models for experiment-1                                        #
#                                                                                                #
##################################################################################################
########################################################################
#                                                                      #
#               2.1 Creating variables of interest                     #
#                                                                      #
########################################################################


### 1 Load the data and assign the unique ID to participants. Generate the aggregated DV for all scenarios.

#creating datset backup
mlm_data_pl <- poland_final_data

#creating datset backup
#assinging ID to each participant
mlm_data_pl$ID <- 1:nrow(mlm_data_pl)
mlm_data_pl$agg_s1 <- (mlm_data_pl$FAIR_S1 + mlm_data_pl$BIAS_S1 + mlm_data_pl$TRUST_S1 + mlm_data_pl$LEGIT_S1)/4
mlm_data_pl$agg_s2 <- (mlm_data_pl$FAIR_S2 + mlm_data_pl$BIAS_S2 + mlm_data_pl$TRUST_S2 + mlm_data_pl$LEGIT_S2)/4
mlm_data_pl$agg_s3 <- (mlm_data_pl$FAIR_S3 + mlm_data_pl$BIAS_S3 + mlm_data_pl$TRUST_S3 + mlm_data_pl$LEGIT_S3)/4
mlm_data_pl$agg_s4 <- (mlm_data_pl$FAIR_S4 + mlm_data_pl$BIAS_S4 + mlm_data_pl$TRUST_S4 + mlm_data_pl$LEGIT_S4)/4


### 2	Convert the wide data format to long data as MLM uses wide data format. Since, we have four scenarios we
###       consider each scenario as level. hence the data has four levels for study1.

mlm_data_pl_L <- melt(mlm_data_pl, id.vars = c("ID"), measure.vars = c("agg_s1", "agg_s2", "agg_s3", "agg_s4"))

# for PL 1st dataset
mlm_data_pl_L2 <- melt(mlm_data_pl, id.vars = c("ID", "CONDITION") , measure.vars = c("agg_s1", "agg_s2", "agg_s3", "agg_s4"))
# for PL second dataset
#mlm_data_pl_L2 <- melt(mlm_data_pl, id.vars = c("ID", "conditionexp1") , measure.vars = c("agg_s1", "agg_s2", "agg_s3", "agg_s4"))


mlm_data_pl_L2$scenario <- recode(mlm_data_pl_L2$variable, '"agg_s1" = 1; "agg_s2" = 2;
"agg_s3" = 3; "agg_s4" = 4;')

#Create dummy variables for conditions and scenarios (S1: moderating discussion, S2:giving civil reminder, S3:presenting facts, S4: presenting news)

mlm_data_pl_L2$AI <- ifelse(mlm_data_pl_L2$CONDITION == 1, 1, 0)
mlm_data_pl_L2$AI_Human <- ifelse(mlm_data_pl_L2$CONDITION == 3, 1, 0)
mlm_data_pl_L2$Human <- ifelse(mlm_data_pl_L2$CONDITION == 2, 1, 0)
mlm_data_pl_L2$S1 <- ifelse(mlm_data_pl_L2$scenario == 1, 1, 0)
mlm_data_pl_L2$S2 <- ifelse(mlm_data_pl_L2$scenario == 2, 1, 0)
mlm_data_pl_L2$S3 <- ifelse(mlm_data_pl_L2$scenario == 3, 1, 0)
mlm_data_pl_L2$S4 <- ifelse(mlm_data_pl_L2$scenario == 4, 1, 0)

########################################################################
#                                                                      #
#            2.2   Multilevel models and main analysis                 #
#                                                                      #
########################################################################


###Multilevel Models for testing H1 and H2

##Human and news recommendation as the reference categories (H1, Appendix 8. Table 1) 

model_e1_pl <- lme(value ~ AI + AI_Human + S1 + S2 + S3 +
                     AI*S1 + AI*S2 + AI*S3 + 
                     AI_Human*S1 + AI_Human*S2 + AI_Human*S3,
                   random =~ 1|ID, data=mlm_data_pl_L2, control = lmeControl(opt = "optim"))
summary(model_e1_pl2)
anova(model_e1_pl)
r.squaredGLMM(model_e1_pl)

##Estimated marginal means (Appendix 8. Table 2)
# The Appendix 8. Table 2 report the the second PL dataset
#Recoding scenario instead of using dummy variables to make visualization easier (coding news as reference cat).
mlm_data_pl_L2$scenario_rec <- ifelse(mlm_data_pl_L2$scenario == 4, 1, 
                                      ifelse(mlm_data_pl_L2$scenario == 1, 4,
                                             ifelse(mlm_data_pl_L2$scenario == 2, 2, 3)))

#Recoding condition instead of using dummy variables to make visualization easier (coding human as reference cat).

mlm_data_pl_L2$condition_rec <- ifelse(mlm_data_pl_L2$CONDITION == 2, 0, 
                                       ifelse(mlm_data_pl_L2$CONDITION == 1, 1, 2))



#Multilevel Model
model_e1_pl <- lme(value ~ factor(condition_rec)*factor(scenario_rec),
                   random =~ 1|ID, data=mlm_data_pl_L2, control = lmeControl(opt = "optim"))


summary(model_e1_pl)
anova(model_e1_pl)

### Getting the estimated marginal means
emm_pl = emmeans::emmeans( model_e1_pl, ~ factor(condition_rec) * factor(scenario_rec))
emm_pl


#AI as the reference category -- RQ1 (Figure 2)

model_e1_pl_AI <- lme(value ~ Human + AI_Human + S1 + S2 + S3 +
                        Human*S1 + Human*S2 + Human*S3 + 
                        AI_Human*S1 + AI_Human*S2 + AI_Human*S3,
                      random =~ 1|ID, data=mlm_data_pl_L2, control = lmeControl(opt = "optim"))
summary(model_e1_pl_AI)
anova(model_e1_pl_AI)

## Reestimated models with present facts as the reference category for Figure 2

###Multilevel Models for Figure 2

#Human as the reference category 

model_e1_pl_facts <- lme(value ~ AI + AI_Human + S1 + S2 + S4 +
                           AI*S1 + AI*S2 + AI*S4 + 
                           AI_Human*S1 + AI_Human*S2 + AI_Human*S4 ,
                         random =~ 1|ID, data=mlm_data_pl_L2)
summary(model_e1_pl_facts)
anova(model_e1_pl_facts)


model_e1_pl_AIfacts <- lme(value ~ Human + AI_Human + S1 + S2 + S4 +
                             Human*S1 + Human*S2 + Human*S4 + 
                             AI_Human*S1 + AI_Human*S2 + AI_Human*S4 ,
                           random =~ 1|ID, data=mlm_data_pl_L2)

summary(model_e1_pl_AIfacts)
anova(model_e1_pl_AIfacts)


########################################################################
#                                                                      #
#     2.3  Multilevel models for pre-registered outcomes               #
#                     (Appendix 10. Experiment 1)                      #
#                                                                      #
########################################################################


####################################### FAIR ################################################################################################
### 3	Convert the wide data format to long data as MLM uses wide data format. Since, we have four scenarios we
###       consider each scenario as level. hence the data has four levels for experiment 1/

mlm_data_pl_L2 <- melt(mlm_data_pl, id.vars = c("ID", "CONDITION") , measure.vars = c("FAIR_S1", "FAIR_S2", "FAIR_S3", "FAIR_S4"))

library(car)
mlm_data_pl_L2$scenario <- recode(mlm_data_pl_L2$variable, '"FAIR_S1" = 1; "FAIR_S2" = 2;
"FAIR_S3" = 3; "FAIR_S4" = 4;')

#Create dummy variables for conditions

mlm_data_pl_L2$AI <- ifelse(mlm_data_pl_L2$CONDITION == 1, 1, 0)
mlm_data_pl_L2$AI_Human <- ifelse(mlm_data_pl_L2$CONDITION == 3, 1, 0)
mlm_data_pl_L2$Human <- ifelse(mlm_data_pl_L2$CONDITION == 2, 1, 0)
mlm_data_pl_L2$S1 <- ifelse(mlm_data_pl_L2$scenario == 1, 1, 0)
mlm_data_pl_L2$S2 <- ifelse(mlm_data_pl_L2$scenario == 2, 1, 0)
mlm_data_pl_L2$S3 <- ifelse(mlm_data_pl_L2$scenario == 3, 1, 0)
mlm_data_pl_L2$S4 <- ifelse(mlm_data_pl_L2$scenario == 4, 1, 0)

#Multilevel Model
#Multilevel Models for figures
# human as a ref

model_e1_pl_fair <- lme(value ~ AI + AI_Human + S1 + S2 + S3 +
                          AI*S1 + AI*S2 + AI*S3 + 
                          AI_Human*S1 + AI_Human*S2 + AI_Human*S3,
                        random =~ 1|ID, data=mlm_data_pl_L2, control = lmeControl(opt = "optim"))
summary(model_e1_pl_fair)
anova(model_e1_pl_fair)

#Recoding scenario instead of using dummy variables to make visualization easier (coding news as reference cat).
mlm_data_pl_L2$scenario_rec <- ifelse(mlm_data_pl_L2$scenario == 4, 1, 
                                      ifelse(mlm_data_pl_L2$scenario == 1, 4,
                                             ifelse(mlm_data_pl_L2$scenario == 2, 2, 3)))

#Recoding condition instead of using dummy variables to make visualization easier (coding human as reference cat).

mlm_data_pl_L2$condition_rec <- ifelse(mlm_data_pl_L2$CONDITION == 2, 0, 
                                       ifelse(mlm_data_pl_L2$CONDITION == 1, 1, 2))

### Condition 0 = Human, Condition 1 = AI, Condition 2 = Mixed
### Scenario 1 = News, Scenario 2 = Civil Reminder, Scenario 3 = Present Facts, Scenario 4 = Moderate


#Multilevel Model
model_e1_pl_mfair <- lme(value ~ factor(condition_rec)*factor(scenario_rec),
                         random =~ 1|ID, data=mlm_data_pl_L2, control = lmeControl(opt = "optim"))
summary(model_e1_pl_mfair)
anova(model_e1_pl_mfair)
vif.lme(model_e1_pl_mfair)

### Getting the estimated marginal means
emm1_fair = emmeans::emmeans(model_e1_pl_mfair, ~ factor(condition_rec) * factor(scenario_rec))
emm1_fair

####################################### BIAS ################################################################################################
### 3	Convert the wide data format to long data as MLM uses wide data format. Since, we have four scenarios we
###       consider each scenario as level. hence the data has four levels for experiment 1

mlm_data_pl_L2 <- melt(mlm_data_pl, id.vars = c("ID", "CONDITION") , measure.vars = c("BIAS_S1", "BIAS_S2", "BIAS_S3", "BIAS_S4"))

mlm_data_pl_L2$scenario <- recode(mlm_data_pl_L2$variable, '"BIAS_S1" = 1; "BIAS_S2" = 2;
"BIAS_S3" = 3; "BIAS_S4" = 4;')
#Create dummy variables for conditions

mlm_data_pl_L2$AI <- ifelse(mlm_data_pl_L2$CONDITION == 1, 1, 0)
mlm_data_pl_L2$AI_Human <- ifelse(mlm_data_pl_L2$CONDITION == 3, 1, 0)
mlm_data_pl_L2$Human <- ifelse(mlm_data_pl_L2$CONDITION == 2, 1, 0)
mlm_data_pl_L2$S1 <- ifelse(mlm_data_pl_L2$scenario == 1, 1, 0)
mlm_data_pl_L2$S2 <- ifelse(mlm_data_pl_L2$scenario == 2, 1, 0)
mlm_data_pl_L2$S3 <- ifelse(mlm_data_pl_L2$scenario == 3, 1, 0)
mlm_data_pl_L2$S4 <- ifelse(mlm_data_pl_L2$scenario == 4, 1, 0)

#Multilevel Model
#Multilevel Models for figures
# human as a ref

model_e1_pl_bias <- lme(value ~ AI + AI_Human + S1 + S2 + S3 +
                          AI*S1 + AI*S2 + AI*S3 + 
                          AI_Human*S1 + AI_Human*S2 + AI_Human*S3,
                        random =~ 1|ID, data=mlm_data_pl_L2, control = lmeControl(opt = "optim"))
summary(model_e1_pl_bias)
anova(model_e1_pl_bias)

#Recoding scenario instead of using dummy variables to make visualization easier (coding news as reference cat).
mlm_data_pl_L2$scenario_rec <- ifelse(mlm_data_pl_L2$scenario == 4, 1, 
                                      ifelse(mlm_data_pl_L2$scenario == 1, 4,
                                             ifelse(mlm_data_pl_L2$scenario == 2, 2, 3)))

#Recoding condition instead of using dummy variables to make visualization easier (coding human as reference cat).

mlm_data_pl_L2$condition_rec <- ifelse(mlm_data_pl_L2$CONDITION == 2, 0, 
                                       ifelse(mlm_data_pl_L2$CONDITION == 1, 1, 2))

### Condition 0 = Human, Condition 1 = AI, Condition 2 = Mixed
### Scenario 1 = News, Scenario 2 = Civil Reminder, Scenario 3 = Present Facts, Scenario 4 = Moderate


#Multilevel Model
model_e1_pl_mbias <- lme(value ~ factor(condition_rec)*factor(scenario_rec),
                         random =~ 1|ID, data=mlm_data_pl_L2, control = lmeControl(opt = "optim"))
summary(model_e1_pl_mbias)
anova(model_e1_pl_mbias)
vif.lme(model_e1_pl_mbias)

### Getting the estimated marginal means
emm1_bias = emmeans::emmeans(model_e1_pl_mbias, ~ factor(condition_rec) * factor(scenario_rec))
emm1_bias

####################################### TRUST ################################################################################################
### 3	Convert the wide data format to long data as MLM uses wide data format. Since, we have four scenarios we
###       consider each scenario as level. hence the data has four levels for experiment 1

mlm_data_pl_L2 <- melt(mlm_data_pl, id.vars = c("ID", "CONDITION") , measure.vars = c("TRUST_S1", "TRUST_S2", "TRUST_S3", "TRUST_S4"))

mlm_data_pl_L2$scenario <- recode(mlm_data_pl_L2$variable, '"TRUST_S1" = 1; "TRUST_S2" = 2;
"TRUST_S3" = 3; "TRUST_S4" = 4;')

#Create dummy variables for conditions
mlm_data_pl_L2$AI <- ifelse(mlm_data_pl_L2$CONDITION == 1, 1, 0)
mlm_data_pl_L2$AI_Human <- ifelse(mlm_data_pl_L2$CONDITION == 3, 1, 0)
mlm_data_pl_L2$Human <- ifelse(mlm_data_pl_L2$CONDITION == 2, 1, 0)
mlm_data_pl_L2$S1 <- ifelse(mlm_data_pl_L2$scenario == 1, 1, 0)
mlm_data_pl_L2$S2 <- ifelse(mlm_data_pl_L2$scenario == 2, 1, 0)
mlm_data_pl_L2$S3 <- ifelse(mlm_data_pl_L2$scenario == 3, 1, 0)
mlm_data_pl_L2$S4 <- ifelse(mlm_data_pl_L2$scenario == 4, 1, 0)

#Multilevel Model
#Multilevel Models for figures
# human as a ref

model_e1_pl_trust <- lme(value ~ AI + AI_Human + S1 + S2 + S3 +
                           AI*S1 + AI*S2 + AI*S3 + 
                           AI_Human*S1 + AI_Human*S2 + AI_Human*S3,
                         random =~ 1|ID, data=mlm_data_pl_L2, control = lmeControl(opt = "optim"))
summary(model_e1_pl_trust)
anova(model_e1_pl_trust)

#Recoding scenario instead of using dummy variables to make visualization easier (coding news as reference cat).
mlm_data_pl_L2$scenario_rec <- ifelse(mlm_data_pl_L2$scenario == 4, 1, 
                                      ifelse(mlm_data_pl_L2$scenario == 1, 4,
                                             ifelse(mlm_data_pl_L2$scenario == 2, 2, 3)))

#Recoding condition instead of using dummy variables to make visualization easier (coding human as reference cat).

mlm_data_pl_L2$condition_rec <- ifelse(mlm_data_pl_L2$CONDITION == 2, 0, 
                                       ifelse(mlm_data_pl_L2$CONDITION == 1, 1, 2))

### Condition 0 = Human, Condition 1 = AI, Condition 2 = Mixed
### Scenario 1 = News, Scenario 2 = Civil Reminder, Scenario 3 = Present Facts, Scenario 4 = Moderate


#Multilevel Model
model_e1_pl_mtrust <- lme(value ~ factor(condition_rec)*factor(scenario_rec),
                          random =~ 1|ID, data=mlm_data_pl_L2, control = lmeControl(opt = "optim"))
summary(model_e1_pl_mtrust)
anova(model_e1_pl_mtrust)
vif.lme(model_e1_pl_mtrust)

### Getting the estimated marginal means
emm1_trust = emmeans::emmeans(model_e1_pl_mtrust, ~ factor(condition_rec) * factor(scenario_rec))
emm1_trust

####################################### AGREE ################################################################################################
mlm_data_pl_L2 <- melt(mlm_data_pl, id.vars = c("ID", "CONDITION") , measure.vars = c("LEGIT_S1", "LEGIT_S2", "LEGIT_S3", "LEGIT_S4"))

mlm_data_pl_L2$scenario <- recode(mlm_data_pl_L2$variable, '"LEGIT_S1" = 1; "LEGIT_S2" = 2;
"LEGIT_S3" = 3; "LEGIT_S4" = 4;')

#Create dummy variables for conditions
mlm_data_pl_L2$AI <- ifelse(mlm_data_pl_L2$CONDITION == 1, 1, 0)
mlm_data_pl_L2$AI_Human <- ifelse(mlm_data_pl_L2$CONDITION == 3, 1, 0)
mlm_data_pl_L2$Human <- ifelse(mlm_data_pl_L2$CONDITION == 2, 1, 0)
mlm_data_pl_L2$S1 <- ifelse(mlm_data_pl_L2$scenario == 1, 1, 0)
mlm_data_pl_L2$S2 <- ifelse(mlm_data_pl_L2$scenario == 2, 1, 0)
mlm_data_pl_L2$S3 <- ifelse(mlm_data_pl_L2$scenario == 3, 1, 0)
mlm_data_pl_L2$S4 <- ifelse(mlm_data_pl_L2$scenario == 4, 1, 0)

#Multilevel Model
#Multilevel Models for figures
# human as a ref

model_e1_pl_agree <- lme(value ~ AI + AI_Human + S1 + S2 + S3 +
                           AI*S1 + AI*S2 + AI*S3 + 
                           AI_Human*S1 + AI_Human*S2 + AI_Human*S3,
                         random =~ 1|ID, data=mlm_data_pl_L2, control = lmeControl(opt = "optim"))
summary(model_e1_pl_agree)
anova(model_e1_pl_agree)

#Recoding scenario instead of using dummy variables to make visualization easier (coding news as reference cat).
mlm_data_pl_L2$scenario_rec <- ifelse(mlm_data_pl_L2$scenario == 4, 1, 
                                      ifelse(mlm_data_pl_L2$scenario == 1, 4,
                                             ifelse(mlm_data_pl_L2$scenario == 2, 2, 3)))

#Recoding condition instead of using dummy variables to make visualization easier (coding human as reference cat).

mlm_data_pl_L2$condition_rec <- ifelse(mlm_data_pl_L2$CONDITION == 2, 0, 
                                       ifelse(mlm_data_pl_L2$CONDITION == 1, 1, 2))

### Condition 0 = Human, Condition 1 = AI, Condition 2 = Mixed
### Scenario 1 = News, Scenario 2 = Civil Reminder, Scenario 3 = Present Facts, Scenario 4 = Moderate


#Multilevel Model
model_e1_pl_magree <- lme(value ~ factor(condition_rec)*factor(scenario_rec),
                          random =~ 1|ID, data=mlm_data_pl_L2, control = lmeControl(opt = "optim"))
summary(model_e1_pl_magree)
anova(model_e1_pl_magree)
vif.lme(model_e1_pl_magree)

### Getting the estimated marginal means
emm1_agree = emmeans::emmeans(model_e1_pl_magree, ~ factor(condition_rec) * factor(scenario_rec))
emm1_agree

########################################################################
#                                                                      #
#  Creating variables of interest                                      #
#                                                                      #
########################################################################

# DV creation code - please note that a different statistical software (SPSS) was initially used to create the DVs 
# code below with # at the start of the lines

#corr agreeecon agreegender agreehealth /missing listwise.
#corr qualityecon qualitygender qualityhealth /missing listwise.
#corr credibilityecon credibilitygender credibilityhealth /missing listwise.
#corr biasecon biasgender biashealth /missing listwise.
#RELIABILITY
#/VARIABLES=credibilityecon credibilitygender credibilityhealth
#/SCALE('ALL VARIABLES') ALL
#/MODEL=ALPHA
#/STATISTICS=SCALE
#/SUMMARY=TOTAL.

#FACTOR
#/VARIABLES credibilityecon credibilitygender credibilityhealth
#/MISSING LISTWISE 
#/ANALYSIS credibilityecon credibilitygender credibilityhealth
#/PRINT INITIAL EXTRACTION ROTATION
#/CRITERIA MINEIGEN(1) ITERATE(25)
/EXTRACTION PC
/CRITERIA ITERATE(25)
/ROTATION VARIMAX
/METHOD=CORRELATION.

compute sp2mean.agree = mean.1 (agreeecon ,agreegender ,agreehealth).
compute sp2mean.quality = mean.1 (qualityecon, qualitygender ,qualityhealth).
compute sp2mean.cred = mean.1 (credibilityecon, credibilitygender ,credibilityhealth ).
compute sp2mean.biased = mean.1 (biasecon ,biasgender, biashealth).




# temporary code to check the model's outcomes.
setwd('C:/Users/arti/Documents/AI-perception/scripts/MLM/AI Paper-Joao/')
mlm_data_pl_temp <- read.csv("Poland_final_data_V2.csv")

#creating dataset backup
mlm_data_pl <- mlm_data_pl_temp
mlm_data_pl$govstr <- scale(abs(mlm_data_pl$govsuppopp-5))
mlm_data_pl$ideostr <- scale(abs(mlm_data_pl$ideology-5))
mlm_data_pl$gender_str = scale(abs(mlm_data_pl$attgender1 + mlm_data_pl$attgender2 + mlm_data_pl$attgender3)/3 - 6.33)
mlm_data_pl$health_str = scale(abs(mlm_data_pl$atthealth1 + mlm_data_pl$atthealth2 + mlm_data_pl$atthealth3)/3 - 6.33)
mlm_data_pl$econ_str = scale(abs(mlm_data_pl$govecon1 + mlm_data_pl$govecon2 +  mlm_data_pl$govecon3)/3 - 6.33)

### Computing party identity
mlm_data_pl$suppop  <- ifelse(mlm_data_pl$govsuppopp < 5, 0,
                              ifelse(mlm_data_pl$govsuppopp == 5, 1, 2))
mlm_data_pl$opponent <- ifelse(mlm_data_pl$suppop == 0, 1, 0)
mlm_data_pl$neutral <- ifelse(mlm_data_pl$suppop == 1, 1, 0)
mlm_data_pl$supporter <- ifelse(mlm_data_pl$suppop == 2, 1, 0)

mlm_data_pl$avidopp <- (mlm_data_pl$identOPP1+mlm_data_pl$identOPP2+mlm_data_pl$identOPP3+mlm_data_pl$identOPP4)/4
mlm_data_pl$avidsupp <- (mlm_data_pl$identSUPP1+mlm_data_pl$identSUPP2+mlm_data_pl$identSUPP3+mlm_data_pl$identSUPP4)/4
mlm_data_pl$avidneu <- (mlm_data_pl$identNEUT1+mlm_data_pl$identNEUT2+mlm_data_pl$identNEUT3+mlm_data_pl$identNEUT4)/4

mlm_data_pl$avidopp[is.na(mlm_data_pl$avidopp)] <- 0
mlm_data_pl$avidsupp[is.na(mlm_data_pl$avidsupp)] <- 0
mlm_data_pl$avidneu[is.na(mlm_data_pl$avidneu)] <- 0

mlm_data_pl$avgppid <- mlm_data_pl$avidopp + mlm_data_pl$avidsupp + mlm_data_pl$avidneu

### 2 Load the data and assign the unique ID to participants. Generate the aggregated DV for all scenarios.

#creating datset backup
#assinging ID to each participant
mlm_data_pl$ID <- 1:nrow(mlm_data_pl)
mlm_data_pl$agg_s1 <- (mlm_data_pl$agreeecon + mlm_data_pl$qualityecon + mlm_data_pl$credibilityecon + mlm_data_pl$biasecon)/4
mlm_data_pl$agg_s2 <- (mlm_data_pl$agreegender + mlm_data_pl$qualitygender + mlm_data_pl$credibilitygender + mlm_data_pl$biasgender)/4
mlm_data_pl$agg_s3 <- (mlm_data_pl$agreehealth + mlm_data_pl$qualityhealth + mlm_data_pl$credibilityhealth + mlm_data_pl$biashealth)/4


### 3	Convert the wide data format to long data as MLM uses wide data format. Since, we have three issues we
###       consider each issue as level. hence the data has three levels for study 2.

#mlm_data_pl_L <- melt(mlm_data_pl, id.vars = c("ID", "conditionexp2", "govstr", "health_str", "gender_str", "econ_str", "avgppid"), measure.vars = c("agg_s1", "agg_s2", "agg_s3"))
mlm_data_pl_L2 <- melt(mlm_data_pl, id.vars = c("ID", "CONDITION", "conditionexp2", "ideostr","govstr", "health_str", "gender_str", "econ_str", "avgppid") , measure.vars = c("agg_s1", "agg_s2", "agg_s3"))

library(car)
mlm_data_pl_L2$scenario <- recode(mlm_data_pl_L2$variable, '"agg_s1" = 1; "agg_s2" = 2;"agg_s3" = 3;')

#Create dummy variables for conditions

mlm_data_pl_L2$AI <- ifelse(mlm_data_pl_L2$CONDITION == 1, 1, 0)
mlm_data_pl_L2$Human <- ifelse(mlm_data_pl_L2$CONDITION == 2, 1, 0)
mlm_data_pl_L2$AI_Human <- ifelse(mlm_data_pl_L2$CONDITION == 3, 1, 0)
mlm_data_pl_L2$econ <- ifelse(mlm_data_pl_L2$scenario == 1, 1, 0)
mlm_data_pl_L2$gender <- ifelse(mlm_data_pl_L2$scenario == 2, 1, 0)
mlm_data_pl_L2$health <- ifelse(mlm_data_pl_L2$scenario == 3, 1, 0)
mlm_data_pl_L2$conditionexp2 <- ifelse(mlm_data_pl_L2$conditionexp2 == 2, 1, 0)


########################################################################
#                                                                      #
#        MLM models main analysis (Dr. Joao kindly check this section) #
#                                                                      #
########################################################################
#Multilevel Models (Appendix 9 table 1)
model_e2_pl_H3 <- lme(value ~ conditionexp2 + health + gender +
                        health*conditionexp2 + gender*conditionexp2,
                      random =~ 1|ID, data=mlm_data_pl_L2, control = lmeControl(opt = "optim"))
summary(model_e2_pl_H3)
anova(model_e2_pl_H3)
vif.lme(model_e2_pl_H3)


#Multilevel Models for figures

#Multilevel baseline Model H4 -- ONLY ON THE COUNTER-ATT CONDITION TO AVOID three-way interactions and high values of 
# multicollinearity VIF = 8.15, although still below 10)

## AI as the reference (Appendix 9 table 2)
model_e2_pl_counter <- lme(value ~ Human + AI_Human + health + gender + +
                             Human*health + Human*gender +
                             AI_Human*health + AI_Human*gender,
                           random =~1|ID, data=mlm_data_pl_L2[mlm_data_pl_L2$conditionexp2== 1,], control = lmeControl(opt = "optim"))
summary(model_e2_pl_counter)
anova(model_e2_pl_counter)
vif.lme(model_e2_pl_counter)
car::vifmodel_e2_pl_counter

## mixed as the reference

model_e2_pl_mix <- lme(value ~ Human + AI + health + gender + +
                         Human*health + Human*gender +
                         AI*health + AI*gender,
                       random =~1|ID, data=mlm_data_pl_L2[mlm_data_pl_L2$conditionexp2== 1,], control = lmeControl(opt = "optim"))
summary(model_e2_pl_mix)
anova(model_e2_pl_mix)
vif.lme(model_e2_pl_mix)
car::vifmodel_e2_pl_mix

## EXPLORATORY ANALYSIS -- to test if people more open to pro/counter from ai/human/mixed
#AI as ref cat

model_e2_pl_test <- lme(value ~ Human + AI_Human + conditionexp2 +
                          Human*conditionexp2 + AI_Human*conditionexp2,
                        random =~ 1|ID, data=mlm_data_pl_L2, control = lmeControl(opt = "optim"))
summary(model_e2_pl_test)
anova(model_e2_pl_test)
vif.lme(model_e2_pl_test)

#mixed as ref cat

model_e2_pl_test2 <- lme(value ~ AI + AI_Human + conditionexp2 +
                           AI*conditionexp2 + AI_Human*conditionexp2,
                         random =~ 1|ID, data=mlm_data_pl_L2, control = lmeControl(opt = "optim"))
summary(model_e2_pl_test2)
anova(model_e2_pl_test2)
vif.lme(model_e2_pl_test2)


# H7 Heterogeneous effects (Appendix 12. Table 1)

model_e2_pl_avgppid <- lme(value ~ Human + AI_Human + health + gender + 
                             Human*health + Human*gender +
                             AI_Human*health + AI_Human*gender +
                             avgppid + avgppid*Human + avgppid*AI_Human + 
                             avgppid*Human*gender + avgppid*AI_Human*gender + avgppid*Human*health + avgppid*AI_Human*health,
                           random =~ 1|ID, data=mlm_data_pl_L2[mlm_data_pl_L2$conditionexp2 == 1,], 
                           control = lmeControl(opt = "optim"))
summary(model_e2_pl_avgppid)
anova(model_e2_pl_avgppid)
vif.lme(model_e2_pl_avgppid)

## Ideology strength as a moderator (Appendix 13 table 1)
model_e2_pl_ideostr <- lme(value ~ Human + AI_Human + health + gender + 
                             Human*health + Human*gender +
                             AI_Human*health + AI_Human*gender +
                             ideostr + ideostr*Human + ideostr*AI_Human +
                             ideostr*Human*gender + ideostr*AI_Human*gender + ideostr*Human*health + ideostr*AI_Human*health,
                           random =~ 1|ID, data=mlm_data_pl_L2[mlm_data_pl_L2$conditionexp2 == 1,], 
                           control = lmeControl(opt = "optim"))
summary(model_e2_pl_ideostr)

## Openness to Counter-Attitudinal Information with support/opposition toward government strength in PL
# (Appendix 13 table 4)
model_e2_pl_govstr <- lme(value ~ Human + AI_Human + health + gender + 
                            Human*health + Human*gender +
                            AI_Human*health + AI_Human*gender +
                            govstr + govstr*Human + govstr*AI_Human +
                            govstr*Human*gender + govstr*AI_Human*gender + govstr*Human*health + govstr*AI_Human*health,
                          random =~ 1|ID, data=mlm_data_pl_L2[mlm_data_pl_L2$conditionexp2 == 1,], 
                          control = lmeControl(opt = "optim"))
summary(model_e2_pl_govstr)
anova(model_e2_pl_govstr)
vif.lme(model_e2_pl_govstr)

# H7 Heterogeneous effects issue attitudes (Interacting scenarios with issue att_str, only for counter condition)
#(Appendix 13 table 2) (Dr. Joao Could you check the if I am using the correct model here as I am getting different output that what we have in appendix)

model_e2_pl_att <- lme(value ~ Human + AI_Human + health + gender + 
                         Human*health + Human*gender +  AI_Human*health + AI_Human*gender +
                         gender_str + health_str + econ_str + Human*econ_str + AI_Human*econ_str +
                         Human*health_str + AI_Human*health_str + Human*gender_str + AI_Human*gender_str,
                       random =~ 1|ID, data=mlm_data_pl_L2[mlm_data_pl_L2$conditionexp2 == 1,], 
                       control = lmeControl(opt = "optim"))
summary(model_e2_pl_att)
anova(model_e2_pl_att)
vif.lme(model_e2_pl_att)


########################################################################
#                                                                      #
#       MLM models per DV for the experiment2 (Appendix 11 table 1)    #
#                                                                      #
########################################################################
# creating table per DV for experiment 2 for Poland
# Creating dataset for Agree
mlm_data_pl_e2_agree <- melt(mlm_data_pl, id.vars = c("ID", "conditionexp1", "conditionexp2") , measure.vars = c("agreeecon", "agreegender", "agreehealth"))
mlm_data_pl_e2_agree$scenario <- recode(mlm_data_pl_e2_agree$variable, '"agreeecon" = 1; "agreegender" = 2; "agreehealth" = 3;')


mlm_data_pl_e2_agree$AI <- ifelse(mlm_data_pl_e2_agree$conditionexp1 == 1, 1, 0)
mlm_data_pl_e2_agree$Human <- ifelse(mlm_data_pl_e2_agree$conditionexp1 == 2, 1, 0)
mlm_data_pl_e2_agree$AI_Human <- ifelse(mlm_data_pl_e2_agree$conditionexp1 == 3, 1, 0)
mlm_data_pl_e2_agree$econ <- ifelse(mlm_data_pl_e2_agree$scenario == 1, 1, 0)
mlm_data_pl_e2_agree$gender <- ifelse(mlm_data_pl_e2_agree$scenario == 2, 1, 0)
mlm_data_pl_e2_agree$health <- ifelse(mlm_data_pl_e2_agree$scenario == 3, 1, 0)
mlm_data_pl_e2_agree$counter_pro <- ifelse(mlm_data_pl_e2_agree$conditionexp2 == 2, 1, 0)

model_e2_pl_H4_agree <- lme(value ~ Human + AI_Human + health + gender + +
                              Human*health + Human*gender +
                              AI_Human*health + AI_Human*gender,
                            random =~1|ID, data=mlm_data_pl_e2_agree[mlm_data_pl_e2_agree$conditionexp2== 1,], control = lmeControl(opt = "optim"))
summary(model_e2_pl_H4_agree)

# Creating dataset for Quality
mlm_data_pl_e2_quality <- melt(mlm_data_pl, id.vars = c("ID", "conditionexp1", "conditionexp2") , measure.vars = c("qualityecon", "qualitygender", "qualityhealth"))
mlm_data_pl_e2_quality$scenario <- recode(mlm_data_pl_e2_quality$variable, '"qualityecon" = 1; "qualitygender" = 2; "qualityhealth" = 3;')

mlm_data_pl_e2_quality$AI <- ifelse(mlm_data_pl_e2_quality$conditionexp1 == 1, 1, 0)
mlm_data_pl_e2_quality$Human <- ifelse(mlm_data_pl_e2_quality$conditionexp1 == 2, 1, 0)
mlm_data_pl_e2_quality$AI_Human <- ifelse(mlm_data_pl_e2_quality$conditionexp1 == 3, 1, 0)
mlm_data_pl_e2_quality$econ <- ifelse(mlm_data_pl_e2_quality$scenario == 1, 1, 0)
mlm_data_pl_e2_quality$gender <- ifelse(mlm_data_pl_e2_quality$scenario == 2, 1, 0)
mlm_data_pl_e2_quality$health <- ifelse(mlm_data_pl_e2_quality$scenario == 3, 1, 0)
mlm_data_pl_e2_quality$counter_pro <- ifelse(mlm_data_pl_e2_quality$conditionexp2 == 2, 1, 0)

model_e2_pl_H4_quality <- lme(value ~ Human + AI_Human + health + gender + +
                                Human*health + Human*gender +
                                AI_Human*health + AI_Human*gender,
                              random =~1|ID, data=mlm_data_pl_e2_quality[mlm_data_pl_e2_quality$conditionexp2== 1,], control = lmeControl(opt = "optim"))
summary(model_e2_pl_H4_quality)

# Creating dataset for Credebility
mlm_data_pl_e2_cred <- melt(mlm_data_pl, id.vars = c("ID", "conditionexp1", "conditionexp2") , measure.vars = c("credibilityecon", "credibilitygender", "credibilityhealth"))
mlm_data_pl_e2_cred$scenario <- recode(mlm_data_pl_e2_cred$variable, '"credibilityecon" = 1; "credibilitygender" = 2; "credibilityhealth" = 3;')

mlm_data_pl_e2_cred$AI <- ifelse(mlm_data_pl_e2_cred$conditionexp1 == 1, 1, 0)
mlm_data_pl_e2_cred$Human <- ifelse(mlm_data_pl_e2_cred$conditionexp1 == 2, 1, 0)
mlm_data_pl_e2_cred$AI_Human <- ifelse(mlm_data_pl_e2_cred$conditionexp1 == 3, 1, 0)
mlm_data_pl_e2_cred$econ <- ifelse(mlm_data_pl_e2_cred$scenario == 1, 1, 0)
mlm_data_pl_e2_cred$gender <- ifelse(mlm_data_pl_e2_cred$scenario == 2, 1, 0)
mlm_data_pl_e2_cred$health <- ifelse(mlm_data_pl_e2_cred$scenario == 3, 1, 0)
mlm_data_pl_e2_cred$counter_pro <- ifelse(mlm_data_pl_e2_cred$conditionexp2 == 2, 1, 0)

model_e2_pl_H4_cred <- lme(value ~ Human + AI_Human + health + gender + +
                             Human*health + Human*gender +
                             AI_Human*health + AI_Human*gender,
                           random =~1|ID, data=mlm_data_pl_e2_cred[mlm_data_pl_e2_cred$conditionexp2== 1,], control = lmeControl(opt = "optim"))
summary(model_e2_pl_H4_cred)


# Creating dataset for Bias

mlm_data_pl_e2_bias <- melt(mlm_data_pl, id.vars = c("ID", "conditionexp1", "conditionexp2") , measure.vars = c("biasecon", "biasgender", "biashealth"))
mlm_data_pl_e2_bias$scenario <- recode(mlm_data_pl_e2_bias$variable, '"biasecon" = 1; "biasgender" = 2; "biashealth" = 3;')

mlm_data_pl_e2_bias$AI <- ifelse(mlm_data_pl_e2_bias$conditionexp1 == 1, 1, 0)
mlm_data_pl_e2_bias$Human <- ifelse(mlm_data_pl_e2_bias$conditionexp1 == 2, 1, 0)
mlm_data_pl_e2_bias$AI_Human <- ifelse(mlm_data_pl_e2_bias$conditionexp1 == 3, 1, 0)
mlm_data_pl_e2_bias$econ <- ifelse(mlm_data_pl_e2_bias$scenario == 1, 1, 0)
mlm_data_pl_e2_bias$gender <- ifelse(mlm_data_pl_e2_bias$scenario == 2, 1, 0)
mlm_data_pl_e2_bias$health <- ifelse(mlm_data_pl_e2_bias$scenario == 3, 1, 0)
mlm_data_pl_e2_bias$counter_pro <- ifelse(mlm_data_pl_e2_bias$conditionexp2 == 2, 1, 0)

model_e2_pl_H4_bias <- lme(value ~ Human + AI_Human + health + gender + +
                             Human*health + Human*gender +
                             AI_Human*health + AI_Human*gender,
                           random =~1|ID, data=mlm_data_pl_e2_bias[mlm_data_pl_e2_bias$conditionexp2== 1,], control = lmeControl(opt = "optim"))
summary(model_e2_pl_H4_bias)


########################################################################
#                                                                      #
#   Composite Marginal means for the experiment2                       #
#                                                                      #
########################################################################
mlm_data_pl$govstr <- scale(abs(mlm_data_pl$govsuppopp-5))
mlm_data_pl$gender_str = scale(abs(mlm_data_pl$attgender1 + mlm_data_pl$attgender2 + mlm_data_pl$attgender3)/3 - 6.33)
mlm_data_pl$health_str = scale(abs(mlm_data_pl$atthealth1 + mlm_data_pl$atthealth2 + mlm_data_pl$atthealth3)/3 - 6.33)
mlm_data_pl$econ_str = scale(abs(mlm_data_pl$govecon1 + mlm_data_pl$govecon2 +  mlm_data_pl$govecon3)/3 - 6.33)

### Computing party identity
mlm_data_pl$suppop  <- ifelse(mlm_data_pl$govsuppopp < 5, 0,
                              ifelse(mlm_data_pl$govsuppopp == 5, 1, 2))
mlm_data_pl$opponent <- ifelse(mlm_data_pl$suppop == 0, 1, 0)
mlm_data_pl$neutral <- ifelse(mlm_data_pl$suppop == 1, 1, 0)
mlm_data_pl$supporter <- ifelse(mlm_data_pl$suppop == 2, 1, 0)

mlm_data_pl$avidopp <- (mlm_data_pl$identOPP1+mlm_data_pl$identOPP2+mlm_data_pl$identOPP3+mlm_data_pl$identOPP4)/4
mlm_data_pl$avidsupp <- (mlm_data_pl$identSUPP1+mlm_data_pl$identSUPP2+mlm_data_pl$identSUPP3+mlm_data_pl$identSUPP4)/4
mlm_data_pl$avidneu <- (mlm_data_pl$identNEUT1+mlm_data_pl$identNEUT2+mlm_data_pl$identNEUT3+mlm_data_pl$identNEUT4)/4

mlm_data_pl$avidopp[is.na(mlm_data_pl$avidopp)] <- 0
mlm_data_pl$avidsupp[is.na(mlm_data_pl$avidsupp)] <- 0
mlm_data_pl$avidneu[is.na(mlm_data_pl$avidneu)] <- 0

mlm_data_pl$avgppid <- mlm_data_pl$avidopp + mlm_data_pl$avidsupp + mlm_data_pl$avidneu

### 2 Load the data and assign the unique ID to participants. Generate the aggregated DV for all scenarios.

#creating datset backup
#assinging ID to each participant
mlm_data_pl$ID <- 1:nrow(mlm_data_pl)
mlm_data_pl$agg_s1 <- (mlm_data_pl$agreeecon + mlm_data_pl$qualityecon + mlm_data_pl$credibilityecon + mlm_data_pl$biasecon)/4
mlm_data_pl$agg_s2 <- (mlm_data_pl$agreegender + mlm_data_pl$qualitygender + mlm_data_pl$credibilitygender + mlm_data_pl$biasgender)/4
mlm_data_pl$agg_s3 <- (mlm_data_pl$agreehealth + mlm_data_pl$qualityhealth + mlm_data_pl$credibilityhealth + mlm_data_pl$biashealth)/4


### 3	Convert the wide data format to long data as MLM uses wide data format. Since, we have three issues we
###       consider each issue as level. hence the data has three levels for study 2.

mlm_data_pl_L <- melt(mlm_data_pl, id.vars = c("ID", "conditionexp2", "govstr", "health_str", "gender_str", "econ_str"), measure.vars = c("agg_s1", "agg_s2", "agg_s3"))
mlm_data_pl_L2 <- melt(mlm_data_pl, id.vars = c("ID", "CONDITION", "conditionexp2", "avgppid","govstr", "health_str", "gender_str", "econ_str") , measure.vars = c("agg_s1", "agg_s2", "agg_s3"))

mlm_data_pl_L2$scenario <- recode(mlm_data_pl_L2$variable, '"agg_s1" = 1; "agg_s2" = 2;"agg_s3" = 3;')

#Recoding scenario instead of using dummy variables to make visualization easier (coding trump economy as reference cat).
mlm_data_pl_L2$scenario_rec <- ifelse(mlm_data_pl_L2$scenario == 3, 1, ifelse(mlm_data_pl_L2$scenario == 1, 3, 2))


#Model for H3, H4 and H5 using factor variables instead of dummy (results are the same, only variable names are less clear in output)
#(Appendix-9)
model_e2_pl_emm <- lme(value ~ factor(CONDITION) * factor(scenario_rec) * factor(conditionexp2),
                       random =~ 1|ID, data=mlm_data_pl_L2, control = lmeControl(opt = "optim"))

### Condition 1 = AI, Condition 2 = Human, Condition 3 = Mixed
### Issue 1 = Healthcare, Issue 2 = Gender equality, Issue 3 = Ecnomic policy
### conditionexp2 1 = pro-attitudinal, conditionexp2 2 = counter-attitudinal

summary(model_e2_pl_emm)
anova(model_e2_pl_emm)
vif.lme(model_e2_pl_emm)
### Getting the estimated marginal means
emm1 = emmeans::emmeans(model_e2_pl_emm, ~ factor(CONDITION) * factor(scenario_rec) * factor(conditionexp2))
emm1


# H7 Heterogeneous effects
#Marginal Means for Openness to Counter-Attitudinal Information; Political identity strength as the moderator 
# Appendix 12 table 2
model_e2_pl_avgppid <- lme(value ~ factor(CONDITION) * factor(scenario_rec) * avgppid,
                           random =~ 1|ID, data=mlm_data_pl_L2[mlm_data_pl_L2$conditionexp2 == 1,], 
                           control = lmeControl(opt = "optim"))
summary(model_e2_pl_avgppid)
anova(model_e2_pl_avgppid)
vif.lme(model_e2_pl_avgppid)

### Getting the estimated marginal means.
#NOTE: Because we have a continupus moderator, we need the argument cov.reduce = range
emm2 = emmeans::emmeans(model_e2_pl_avgppid, ~ factor(CONDITION) * factor(scenario_rec) * avgppid, cov.reduce = range)
emm2

#############################################################################################
#                                                                                           #
#                   3  Multilevel models for experiment-2                                   #
#                                                                                           #
#                                                                                           #
#############################################################################################
########################################################################
#                                                                      #
#        3.1  Creating variables of interest                           #
#                                                                      #
########################################################################

# DVs creation code goes here (DVs were created in a different statistical software - SPSS - the code for DV creation is pasted below 
# indicated with #s

#corr agreeecon agreegender agreehealth /missing listwise.
#corr qualityecon qualitygender qualityhealth /missing listwise.
#corr credibilityecon credibilitygender credibilityhealth /missing listwise.
#corr biasecon biasgender biashealth /missing listwise.
#RELIABILITY
#/VARIABLES=credibilityecon credibilitygender credibilityhealth
#/SCALE('ALL VARIABLES') ALL
#/MODEL=ALPHA
#/STATISTICS=SCALE
#/SUMMARY=TOTAL.

#FACTOR
#/VARIABLES credibilityecon credibilitygender credibilityhealth
#/MISSING LISTWISE 
#/ANALYSIS credibilityecon credibilitygender credibilityhealth
#/PRINT INITIAL EXTRACTION ROTATION
#/CRITERIA MINEIGEN(1) ITERATE(25)
#/EXTRACTION PC
#/CRITERIA ITERATE(25)
#/ROTATION VARIMAX
#/METHOD=CORRELATION.

#compute sp2mean.agree = mean.1 (agreeecon ,agreegender ,agreehealth).
#compute sp2mean.quality = mean.1 (qualityecon, qualitygender ,qualityhealth).
#compute sp2mean.cred = mean.1 (credibilityecon, credibilitygender ,credibilityhealth ).
#compute sp2mean.biased = mean.1 (biasecon ,biasgender, biashealth).



# temporary code using previously saved cleaned dataset.
setwd('C:/Users/arti/Documents/AI-perception/scripts/MLM/AI Paper-Joao/')
mlm_data_pl_temp <- read.csv("Poland_final_data_V2.csv")

#creating dataset backup
mlm_data_pl <- mlm_data_pl_temp
mlm_data_pl$govstr <- scale(abs(mlm_data_pl$govsuppopp-5))
mlm_data_pl$ideostr <- scale(abs(mlm_data_pl$ideology-5))
mlm_data_pl$gender_str = scale(abs(mlm_data_pl$attgender1 + mlm_data_pl$attgender2 + mlm_data_pl$attgender3)/3 - 6.33)
mlm_data_pl$health_str = scale(abs(mlm_data_pl$atthealth1 + mlm_data_pl$atthealth2 + mlm_data_pl$atthealth3)/3 - 6.33)
mlm_data_pl$econ_str = scale(abs(mlm_data_pl$govecon1 + mlm_data_pl$govecon2 +  mlm_data_pl$govecon3)/3 - 6.33)

### Computing party identity
mlm_data_pl$suppop  <- ifelse(mlm_data_pl$govsuppopp < 5, 0,
                              ifelse(mlm_data_pl$govsuppopp == 5, 1, 2))
mlm_data_pl$opponent <- ifelse(mlm_data_pl$suppop == 0, 1, 0)
mlm_data_pl$neutral <- ifelse(mlm_data_pl$suppop == 1, 1, 0)
mlm_data_pl$supporter <- ifelse(mlm_data_pl$suppop == 2, 1, 0)

mlm_data_pl$avidopp <- (mlm_data_pl$identOPP1+mlm_data_pl$identOPP2+mlm_data_pl$identOPP3+mlm_data_pl$identOPP4)/4
mlm_data_pl$avidsupp <- (mlm_data_pl$identSUPP1+mlm_data_pl$identSUPP2+mlm_data_pl$identSUPP3+mlm_data_pl$identSUPP4)/4
mlm_data_pl$avidneu <- (mlm_data_pl$identNEUT1+mlm_data_pl$identNEUT2+mlm_data_pl$identNEUT3+mlm_data_pl$identNEUT4)/4

mlm_data_pl$avidopp[is.na(mlm_data_pl$avidopp)] <- 0
mlm_data_pl$avidsupp[is.na(mlm_data_pl$avidsupp)] <- 0
mlm_data_pl$avidneu[is.na(mlm_data_pl$avidneu)] <- 0

mlm_data_pl$avgppid <- mlm_data_pl$avidopp + mlm_data_pl$avidsupp + mlm_data_pl$avidneu

### 2 Load the data and assign the unique ID to participants. Generate the aggregated DV for all scenarios.

#creating datset backup
#assinging ID to each participant
mlm_data_pl$ID <- 1:nrow(mlm_data_pl)
mlm_data_pl$agg_s1 <- (mlm_data_pl$agreeecon + mlm_data_pl$qualityecon + mlm_data_pl$credibilityecon + mlm_data_pl$biasecon)/4
mlm_data_pl$agg_s2 <- (mlm_data_pl$agreegender + mlm_data_pl$qualitygender + mlm_data_pl$credibilitygender + mlm_data_pl$biasgender)/4
mlm_data_pl$agg_s3 <- (mlm_data_pl$agreehealth + mlm_data_pl$qualityhealth + mlm_data_pl$credibilityhealth + mlm_data_pl$biashealth)/4


### 3	Convert the wide data format to long data as MLM uses wide data format. Since, we have three issues we
###       consider each issue as level. hence the data has three levels for study 2.

#mlm_data_pl_L <- melt(mlm_data_pl, id.vars = c("ID", "conditionexp2", "govstr", "health_str", "gender_str", "econ_str", "avgppid"), measure.vars = c("agg_s1", "agg_s2", "agg_s3"))
mlm_data_pl_L2 <- melt(mlm_data_pl, id.vars = c("ID", "CONDITION", "conditionexp2", "ideostr","govstr", "health_str", "gender_str", "econ_str", "avgppid") , measure.vars = c("agg_s1", "agg_s2", "agg_s3"))

library(car)
mlm_data_pl_L2$scenario <- recode(mlm_data_pl_L2$variable, '"agg_s1" = 1; "agg_s2" = 2;"agg_s3" = 3;')

#Create dummy variables for conditions

mlm_data_pl_L2$AI <- ifelse(mlm_data_pl_L2$CONDITION == 1, 1, 0)
mlm_data_pl_L2$Human <- ifelse(mlm_data_pl_L2$CONDITION == 2, 1, 0)
mlm_data_pl_L2$AI_Human <- ifelse(mlm_data_pl_L2$CONDITION == 3, 1, 0)
mlm_data_pl_L2$econ <- ifelse(mlm_data_pl_L2$scenario == 1, 1, 0)
mlm_data_pl_L2$gender <- ifelse(mlm_data_pl_L2$scenario == 2, 1, 0)
mlm_data_pl_L2$health <- ifelse(mlm_data_pl_L2$scenario == 3, 1, 0)
mlm_data_pl_L2$conditionexp2 <- ifelse(mlm_data_pl_L2$conditionexp2 == 2, 1, 0)


########################################################################
#                                                                      #
#        3.2   MLM models main analysis (By Dr. Joao)                  #
#                                                                      #
########################################################################
#Multilevel Models (Appendix 9 table 1)
model_e2_pl_H3 <- lme(value ~ conditionexp2 + health + gender +
                        health*conditionexp2 + gender*conditionexp2,
                      random =~ 1|ID, data=mlm_data_pl_L2, control = lmeControl(opt = "optim"))
summary(model_e2_pl_H3)
anova(model_e2_pl_H3)
vif.lme(model_e2_pl_H3)


#Multilevel Models for figures

#Multilevel baseline Model H4 -- ONLY ON THE COUNTER-ATT CONDITION TO AVOID three-way interactions and high values of 
# multicollinearity VIF = 8.15, although still below 10)

## AI as the reference (Appendix 9 table 2)
model_e2_pl_counter <- lme(value ~ Human + AI_Human + health + gender + +
                             Human*health + Human*gender +
                             AI_Human*health + AI_Human*gender,
                           random =~1|ID, data=mlm_data_pl_L2[mlm_data_pl_L2$conditionexp2== 1,], control = lmeControl(opt = "optim"))
summary(model_e2_pl_counter)
anova(model_e2_pl_counter)
vif.lme(model_e2_pl_counter)
car::vifmodel_e2_pl_counter

## mixed as the reference

model_e2_pl_mix <- lme(value ~ Human + AI + health + gender + +
                         Human*health + Human*gender +
                         AI*health + AI*gender,
                       random =~1|ID, data=mlm_data_pl_L2[mlm_data_pl_L2$conditionexp2== 1,], control = lmeControl(opt = "optim"))
summary(model_e2_pl_mix)
anova(model_e2_pl_mix)
vif.lme(model_e2_pl_mix)
car::vifmodel_e2_pl_mix

## EXPLORATORY ANALYSIS -- to test if people more open to pro/counter from ai/human/mixed
#AI as ref cat

model_e2_pl_test <- lme(value ~ Human + AI_Human + conditionexp2 +
                          Human*conditionexp2 + AI_Human*conditionexp2,
                        random =~ 1|ID, data=mlm_data_pl_L2, control = lmeControl(opt = "optim"))
summary(model_e2_pl_test)
anova(model_e2_pl_test)
vif.lme(model_e2_pl_test)

#mixed as ref cat

model_e2_pl_test2 <- lme(value ~ AI + AI_Human + conditionexp2 +
                           AI*conditionexp2 + AI_Human*conditionexp2,
                         random =~ 1|ID, data=mlm_data_pl_L2, control = lmeControl(opt = "optim"))
summary(model_e2_pl_test2)
anova(model_e2_pl_test2)
vif.lme(model_e2_pl_test2)


# H7 Heterogeneous effects (Appendix 12. Table 1)

model_e2_pl_avgppid <- lme(value ~ Human + AI_Human + health + gender + 
                             Human*health + Human*gender +
                             AI_Human*health + AI_Human*gender +
                             avgppid + avgppid*Human + avgppid*AI_Human + 
                             avgppid*Human*gender + avgppid*AI_Human*gender + avgppid*Human*health + avgppid*AI_Human*health,
                           random =~ 1|ID, data=mlm_data_pl_L2[mlm_data_pl_L2$conditionexp2 == 1,], 
                           control = lmeControl(opt = "optim"))
summary(model_e2_pl_avgppid)
anova(model_e2_pl_avgppid)
vif.lme(model_e2_pl_avgppid)

## Ideology strength as a moderator (Appendix 13 table 1)
model_e2_pl_ideostr <- lme(value ~ Human + AI_Human + health + gender + 
                             Human*health + Human*gender +
                             AI_Human*health + AI_Human*gender +
                             ideostr + ideostr*Human + ideostr*AI_Human +
                             ideostr*Human*gender + ideostr*AI_Human*gender + ideostr*Human*health + ideostr*AI_Human*health,
                           random =~ 1|ID, data=mlm_data_pl_L2[mlm_data_pl_L2$conditionexp2 == 1,], 
                           control = lmeControl(opt = "optim"))
summary(model_e2_pl_ideostr)

## Openness to Counter-Attitudinal Information with support/opposition toward government strength in PL
# (Appendix 13 table 4)
model_e2_pl_govstr <- lme(value ~ Human + AI_Human + health + gender + 
                            Human*health + Human*gender +
                            AI_Human*health + AI_Human*gender +
                            govstr + govstr*Human + govstr*AI_Human +
                            govstr*Human*gender + govstr*AI_Human*gender + govstr*Human*health + govstr*AI_Human*health,
                          random =~ 1|ID, data=mlm_data_pl_L2[mlm_data_pl_L2$conditionexp2 == 1,], 
                          control = lmeControl(opt = "optim"))
summary(model_e2_pl_govstr)
anova(model_e2_pl_govstr)
vif.lme(model_e2_pl_govstr)

# H7 Heterogeneous effects issue attitudes (Interacting scenarios with issue att_str, only for counter condition)
#(Appendix 13 table 2) 

model_e2_pl_att <- lme(value ~ Human + AI_Human + health + gender + 
                         Human*health + Human*gender +  AI_Human*health + AI_Human*gender +
                         gender_str + health_str + econ_str + Human*econ_str + AI_Human*econ_str +
                         Human*health_str + AI_Human*health_str + Human*gender_str + AI_Human*gender_str,
                       random =~ 1|ID, data=mlm_data_pl_L2[mlm_data_pl_L2$conditionexp2 == 1,], 
                       control = lmeControl(opt = "optim"))
summary(model_e2_pl_att)
anova(model_e2_pl_att)
vif.lme(model_e2_pl_att)


########################################################################
#                                                                      #
#   3.3 MLM models per DV for the experiment2 (Appendix 11 table 1)    #
#                                                                      #
########################################################################
# creating table per DV for experiment 2 for Poland
# Creating dataset for Agree
mlm_data_pl_e2_agree <- melt(mlm_data_pl, id.vars = c("ID", "conditionexp1", "conditionexp2") , measure.vars = c("agreeecon", "agreegender", "agreehealth"))
mlm_data_pl_e2_agree$scenario <- recode(mlm_data_pl_e2_agree$variable, '"agreeecon" = 1; "agreegender" = 2; "agreehealth" = 3;')


mlm_data_pl_e2_agree$AI <- ifelse(mlm_data_pl_e2_agree$conditionexp1 == 1, 1, 0)
mlm_data_pl_e2_agree$Human <- ifelse(mlm_data_pl_e2_agree$conditionexp1 == 2, 1, 0)
mlm_data_pl_e2_agree$AI_Human <- ifelse(mlm_data_pl_e2_agree$conditionexp1 == 3, 1, 0)
mlm_data_pl_e2_agree$econ <- ifelse(mlm_data_pl_e2_agree$scenario == 1, 1, 0)
mlm_data_pl_e2_agree$gender <- ifelse(mlm_data_pl_e2_agree$scenario == 2, 1, 0)
mlm_data_pl_e2_agree$health <- ifelse(mlm_data_pl_e2_agree$scenario == 3, 1, 0)
mlm_data_pl_e2_agree$counter_pro <- ifelse(mlm_data_pl_e2_agree$conditionexp2 == 2, 1, 0)

model_e2_pl_H4_agree <- lme(value ~ Human + AI_Human + health + gender + +
                              Human*health + Human*gender +
                              AI_Human*health + AI_Human*gender,
                            random =~1|ID, data=mlm_data_pl_e2_agree[mlm_data_pl_e2_agree$conditionexp2== 1,], control = lmeControl(opt = "optim"))
summary(model_e2_pl_H4_agree)

# Creating dataset for Quality
mlm_data_pl_e2_quality <- melt(mlm_data_pl, id.vars = c("ID", "conditionexp1", "conditionexp2") , measure.vars = c("qualityecon", "qualitygender", "qualityhealth"))
mlm_data_pl_e2_quality$scenario <- recode(mlm_data_pl_e2_quality$variable, '"qualityecon" = 1; "qualitygender" = 2; "qualityhealth" = 3;')

mlm_data_pl_e2_quality$AI <- ifelse(mlm_data_pl_e2_quality$conditionexp1 == 1, 1, 0)
mlm_data_pl_e2_quality$Human <- ifelse(mlm_data_pl_e2_quality$conditionexp1 == 2, 1, 0)
mlm_data_pl_e2_quality$AI_Human <- ifelse(mlm_data_pl_e2_quality$conditionexp1 == 3, 1, 0)
mlm_data_pl_e2_quality$econ <- ifelse(mlm_data_pl_e2_quality$scenario == 1, 1, 0)
mlm_data_pl_e2_quality$gender <- ifelse(mlm_data_pl_e2_quality$scenario == 2, 1, 0)
mlm_data_pl_e2_quality$health <- ifelse(mlm_data_pl_e2_quality$scenario == 3, 1, 0)
mlm_data_pl_e2_quality$counter_pro <- ifelse(mlm_data_pl_e2_quality$conditionexp2 == 2, 1, 0)

model_e2_pl_H4_quality <- lme(value ~ Human + AI_Human + health + gender + +
                                Human*health + Human*gender +
                                AI_Human*health + AI_Human*gender,
                              random =~1|ID, data=mlm_data_pl_e2_quality[mlm_data_pl_e2_quality$conditionexp2== 1,], control = lmeControl(opt = "optim"))
summary(model_e2_pl_H4_quality)

# Creating dataset for Credebility
mlm_data_pl_e2_cred <- melt(mlm_data_pl, id.vars = c("ID", "conditionexp1", "conditionexp2") , measure.vars = c("credibilityecon", "credibilitygender", "credibilityhealth"))
mlm_data_pl_e2_cred$scenario <- recode(mlm_data_pl_e2_cred$variable, '"credibilityecon" = 1; "credibilitygender" = 2; "credibilityhealth" = 3;')

mlm_data_pl_e2_cred$AI <- ifelse(mlm_data_pl_e2_cred$conditionexp1 == 1, 1, 0)
mlm_data_pl_e2_cred$Human <- ifelse(mlm_data_pl_e2_cred$conditionexp1 == 2, 1, 0)
mlm_data_pl_e2_cred$AI_Human <- ifelse(mlm_data_pl_e2_cred$conditionexp1 == 3, 1, 0)
mlm_data_pl_e2_cred$econ <- ifelse(mlm_data_pl_e2_cred$scenario == 1, 1, 0)
mlm_data_pl_e2_cred$gender <- ifelse(mlm_data_pl_e2_cred$scenario == 2, 1, 0)
mlm_data_pl_e2_cred$health <- ifelse(mlm_data_pl_e2_cred$scenario == 3, 1, 0)
mlm_data_pl_e2_cred$counter_pro <- ifelse(mlm_data_pl_e2_cred$conditionexp2 == 2, 1, 0)

model_e2_pl_H4_cred <- lme(value ~ Human + AI_Human + health + gender + +
                             Human*health + Human*gender +
                             AI_Human*health + AI_Human*gender,
                           random =~1|ID, data=mlm_data_pl_e2_cred[mlm_data_pl_e2_cred$conditionexp2== 1,], control = lmeControl(opt = "optim"))
summary(model_e2_pl_H4_cred)


# Creating dataset for Bias

mlm_data_pl_e2_bias <- melt(mlm_data_pl, id.vars = c("ID", "conditionexp1", "conditionexp2") , measure.vars = c("biasecon", "biasgender", "biashealth"))
mlm_data_pl_e2_bias$scenario <- recode(mlm_data_pl_e2_bias$variable, '"biasecon" = 1; "biasgender" = 2; "biashealth" = 3;')

mlm_data_pl_e2_bias$AI <- ifelse(mlm_data_pl_e2_bias$conditionexp1 == 1, 1, 0)
mlm_data_pl_e2_bias$Human <- ifelse(mlm_data_pl_e2_bias$conditionexp1 == 2, 1, 0)
mlm_data_pl_e2_bias$AI_Human <- ifelse(mlm_data_pl_e2_bias$conditionexp1 == 3, 1, 0)
mlm_data_pl_e2_bias$econ <- ifelse(mlm_data_pl_e2_bias$scenario == 1, 1, 0)
mlm_data_pl_e2_bias$gender <- ifelse(mlm_data_pl_e2_bias$scenario == 2, 1, 0)
mlm_data_pl_e2_bias$health <- ifelse(mlm_data_pl_e2_bias$scenario == 3, 1, 0)
mlm_data_pl_e2_bias$counter_pro <- ifelse(mlm_data_pl_e2_bias$conditionexp2 == 2, 1, 0)

model_e2_pl_H4_bias <- lme(value ~ Human + AI_Human + health + gender + +
                             Human*health + Human*gender +
                             AI_Human*health + AI_Human*gender,
                           random =~1|ID, data=mlm_data_pl_e2_bias[mlm_data_pl_e2_bias$conditionexp2== 1,], control = lmeControl(opt = "optim"))
summary(model_e2_pl_H4_bias)


########################################################################
#                                                                      #
#   3.4 Composite Marginal means for the experiment2                   #
#                                                                      #
########################################################################
mlm_data_pl$govstr <- scale(abs(mlm_data_pl$govsuppopp-5))
mlm_data_pl$gender_str = scale(abs(mlm_data_pl$attgender1 + mlm_data_pl$attgender2 + mlm_data_pl$attgender3)/3 - 6.33)
mlm_data_pl$health_str = scale(abs(mlm_data_pl$atthealth1 + mlm_data_pl$atthealth2 + mlm_data_pl$atthealth3)/3 - 6.33)
mlm_data_pl$econ_str = scale(abs(mlm_data_pl$govecon1 + mlm_data_pl$govecon2 +  mlm_data_pl$govecon3)/3 - 6.33)

### Computing party identity
mlm_data_pl$suppop  <- ifelse(mlm_data_pl$govsuppopp < 5, 0,
                              ifelse(mlm_data_pl$govsuppopp == 5, 1, 2))
mlm_data_pl$opponent <- ifelse(mlm_data_pl$suppop == 0, 1, 0)
mlm_data_pl$neutral <- ifelse(mlm_data_pl$suppop == 1, 1, 0)
mlm_data_pl$supporter <- ifelse(mlm_data_pl$suppop == 2, 1, 0)

mlm_data_pl$avidopp <- (mlm_data_pl$identOPP1+mlm_data_pl$identOPP2+mlm_data_pl$identOPP3+mlm_data_pl$identOPP4)/4
mlm_data_pl$avidsupp <- (mlm_data_pl$identSUPP1+mlm_data_pl$identSUPP2+mlm_data_pl$identSUPP3+mlm_data_pl$identSUPP4)/4
mlm_data_pl$avidneu <- (mlm_data_pl$identNEUT1+mlm_data_pl$identNEUT2+mlm_data_pl$identNEUT3+mlm_data_pl$identNEUT4)/4

mlm_data_pl$avidopp[is.na(mlm_data_pl$avidopp)] <- 0
mlm_data_pl$avidsupp[is.na(mlm_data_pl$avidsupp)] <- 0
mlm_data_pl$avidneu[is.na(mlm_data_pl$avidneu)] <- 0

mlm_data_pl$avgppid <- mlm_data_pl$avidopp + mlm_data_pl$avidsupp + mlm_data_pl$avidneu

### 2 Load the data and assign the unique ID to participants. Generate the aggregated DV for all scenarios.

#creating datset backup
#assinging ID to each participant
mlm_data_pl$ID <- 1:nrow(mlm_data_pl)
mlm_data_pl$agg_s1 <- (mlm_data_pl$agreeecon + mlm_data_pl$qualityecon + mlm_data_pl$credibilityecon + mlm_data_pl$biasecon)/4
mlm_data_pl$agg_s2 <- (mlm_data_pl$agreegender + mlm_data_pl$qualitygender + mlm_data_pl$credibilitygender + mlm_data_pl$biasgender)/4
mlm_data_pl$agg_s3 <- (mlm_data_pl$agreehealth + mlm_data_pl$qualityhealth + mlm_data_pl$credibilityhealth + mlm_data_pl$biashealth)/4


### 3	Convert the wide data format to long data as MLM uses wide data format. Since, we have three issues we
###       consider each issue as level. hence the data has three levels for study 2.

mlm_data_pl_L <- melt(mlm_data_pl, id.vars = c("ID", "conditionexp2", "govstr", "health_str", "gender_str", "econ_str"), measure.vars = c("agg_s1", "agg_s2", "agg_s3"))
mlm_data_pl_L2 <- melt(mlm_data_pl, id.vars = c("ID", "CONDITION", "conditionexp2", "avgppid","govstr", "health_str", "gender_str", "econ_str") , measure.vars = c("agg_s1", "agg_s2", "agg_s3"))

mlm_data_pl_L2$scenario <- recode(mlm_data_pl_L2$variable, '"agg_s1" = 1; "agg_s2" = 2;"agg_s3" = 3;')

#Recoding scenario instead of using dummy variables to make visualization easier (coding trump economy as reference cat).
mlm_data_pl_L2$scenario_rec <- ifelse(mlm_data_pl_L2$scenario == 3, 1, ifelse(mlm_data_pl_L2$scenario == 1, 3, 2))


#Model for H3, H4 and H5 using factor variables instead of dummy (results are the same, only variable names are less clear in output)
#(Appendix-9)
model_e2_pl_emm <- lme(value ~ factor(CONDITION) * factor(scenario_rec) * factor(conditionexp2),
                       random =~ 1|ID, data=mlm_data_pl_L2, control = lmeControl(opt = "optim"))

### Condition 1 = AI, Condition 2 = Human, Condition 3 = Mixed
### Issue 1 = Healthcare, Issue 2 = Gender equality, Issue 3 = Ecnomic policy
### conditionexp2 1 = pro-attitudinal, conditionexp2 2 = counter-attitudinal

summary(model_e2_pl_emm)
anova(model_e2_pl_emm)
vif.lme(model_e2_pl_emm)
### Getting the estimated marginal means
emm1 = emmeans::emmeans(model_e2_pl_emm, ~ factor(CONDITION) * factor(scenario_rec) * factor(conditionexp2))
emm1


# H7 Heterogeneous effects
#Marginal Means for Openness to Counter-Attitudinal Information; Political identity strength as the moderator 
# Appendix 12 table 2
model_e2_pl_avgppid <- lme(value ~ factor(CONDITION) * factor(scenario_rec) * avgppid,
                           random =~ 1|ID, data=mlm_data_pl_L2[mlm_data_pl_L2$conditionexp2 == 1,], 
                           control = lmeControl(opt = "optim"))
summary(model_e2_pl_avgppid)
anova(model_e2_pl_avgppid)
vif.lme(model_e2_pl_avgppid)

### Getting the estimated marginal means.
#NOTE: Because we have a continupus moderator, we need the argument cov.reduce = range
emm2 = emmeans::emmeans(model_e2_pl_avgppid, ~ factor(CONDITION) * factor(scenario_rec) * avgppid, cov.reduce = range)
emm2


